{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ppolicar/nfs/miniconda/envs/tsne/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n",
      "/home/ppolicar/nfs/miniconda/envs/tsne/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "import gzip\n",
    "\n",
    "import utils\n",
    "\n",
    "import loompy\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10min 4s, sys: 11.2 s, total: 10min 15s\n",
      "Wall time: 10min 15s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "data = pd.read_table(\"data/GSE63472_P14Retina_merged_digital_expression.txt.gz\", index_col=0)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 48 ms, sys: 8 ms, total: 56 ms\n",
      "Wall time: 58.3 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "cluster_ids = pd.read_table(\"data/retina_clusteridentities.txt\", header=None, index_col=0, squeeze=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0\n",
       "r1_GGCCGCAGTCCG     2\n",
       "r1_CTTGTGCGGGAA     2\n",
       "r1_GCGCAACTGCTC     2\n",
       "r1_GATTGGGAGGCA     2\n",
       "r1_GTGCCGCCTCTC    25\n",
       "Name: 1, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cluster_ids.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reorder\n",
    "cluster_ids = cluster_ids[data.columns.values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(24658, 49300)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only use cells where metadata is available\n",
    "ind = data.columns.isin(cluster_ids.index)\n",
    "data = data.loc[:, ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((24658, 49300), (49300,))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape, cluster_ids.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ~cluster_ids.isna()\n",
    "data = data.loc[:, mask.values]\n",
    "cluster_ids = cluster_ids[mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert not cluster_ids.isna().any(), \"Did not properly remove cells with NaN label\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((24658, 44808), (44808,))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape, cluster_ids.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12.4 s, sys: 0 ns, total: 12.4 s\n",
      "Wall time: 12.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "counts = sp.csr_matrix(data.values)\n",
    "counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2.34 s, sys: 52 ms, total: 2.39 s\n",
      "Wall time: 1.12 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "cpm_counts = utils.calculate_cpm(counts, axis=0)\n",
    "log_counts = utils.log_normalize(cpm_counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Rods                      29400\n",
       "Bipolar cells              6285\n",
       "Amacrine cells             4426\n",
       "Cones                      1868\n",
       "Muller glia                1624\n",
       "Retinal ganglion cells      432\n",
       "Horizontal cells            252\n",
       "Vascular endothelium        252\n",
       "Fibroblasts                  85\n",
       "Microglia                    67\n",
       "Pericytes                    63\n",
       "Astrocytes                   54\n",
       "Name: 1, dtype: int64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cell_types = cluster_ids.astype(object)\n",
    "\n",
    "cell_types.loc[cell_types == 1] = \"Horizontal cells\"\n",
    "cell_types.loc[cell_types == 2] = \"Retinal ganglion cells\"\n",
    "cell_types.loc[cell_types.isin(range(3, 24))] = \"Amacrine cells\"\n",
    "cell_types.loc[cell_types == 24] = \"Rods\"\n",
    "cell_types.loc[cell_types == 25] = \"Cones\"\n",
    "cell_types.loc[cell_types.isin(range(26, 34))] = \"Bipolar cells\"\n",
    "cell_types.loc[cell_types == 34] = \"Muller glia\"\n",
    "cell_types.loc[cell_types == 35] = \"Astrocytes\"\n",
    "cell_types.loc[cell_types == 36] = \"Fibroblasts\"\n",
    "cell_types.loc[cell_types == 37] = \"Vascular endothelium\"\n",
    "cell_types.loc[cell_types == 38] = \"Pericytes\"\n",
    "cell_types.loc[cell_types == 39] = \"Microglia\"\n",
    "\n",
    "cell_types.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dropout based feature selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chosen offset: 0.19\n",
      "CPU times: user 3.87 s, sys: 148 ms, total: 4.02 s\n",
      "Wall time: 1.13 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAD0CAYAAADOibL4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXl4U2X2xz9vkiahTSmFUiyLgFBAKCAIFGQRAUVcEFcQF0AdlxFxGZlxG0XHAXWcnw7ouKGijgijjsqoiAriwiAUlaUoOy1LWdoC3ds0yfv7I7nhJk3b2zYpXd7P8+RpcnOXk4j35Jzzfc8RUkoUCoVCoWhomE61AQqFQqFQhEI5KIVCoVA0SJSDUigUCkWDRDkohUKhUDRIlINSKBQKRYNEOSiFQqFQNEgi5qCEEG8IIY4KIdIreV8IIeYLIXYJITYLIQZGyhaFQqFQND4iGUEtAi6s4v0JQLLvcSvwUgRtUSgUCkUjI2IOSkr5HXCsil0uA96WXn4EWgkhkiJlj0KhUCgaF5ZTeO0OwH7d6wO+bYeCdxRC3Io3ykJE2c+OatOx7leXgIAodzkWj6fu52vgCLwfOVzYLSacbu/31sZhAyCnoAwJxNostIuzA5BXUs7xIicWk6DU5Ql4T8+RvFIKylwIICHWRlyLqAr7aOdqEWWm0OlCSio9X2XHxsdYKXW6KShzhbQzPsZa4dpVvVeT6xo5NtT+dbl+XWxRKAyjdSTS7qVBHYp+2rkzR0rZtqanPZUOyjBSyleBVwFsSckyadrzYTu31eWke84+pRapAwkOKxLILXQCYBKwbOYIAKa9ud6/XSPGamb68C5MSDkZME984Qc8vn/TAoj27QOwJG0/UwZ3YknafnILnZgE/n0BPr3Le63l6YdYtCYDoML5NTsSHFaOFTnxyNB2JjisLJoxJMDeqt6rDv2xk32fYcrgTgG2VXetulw/XJ9Doag12dmICRMya3PoqbwvHwQ66V539G2rN0weD06LlRMtWtbnZRsFZlH19p6JDmKsZhxWM5MHd2LK4E7YzF7nMio5wb//lMGdSHBY6Zno8G8rcrpZmrY/4LyjkhPQLil9+yxak8FLq3eTW+hkadp++rZviUlAclsHoViStp8ipzvk+TU7Jg/uxKjkBEwitJ2TB3cKPm2V71WH/ljNwQbbVt216nL96s6tUDRkRCSbxQohugCfSilTQrx3MTATuAhIBeZLKav9WRfuCArA5ouiKrknNwjMAtwG/lPpo4sYq5nBXeJZvSMH8DqPHokOth8t9L+OtpopcrpJcFhJad+S73bmMCo5gfvH9wr7Z9AiHAEM6hLPlqz8CtHE8vRDvLUmw5+OLHK6MQm4Y3Q3/w0+2NY+HeJY4nNgGzKOI6kYQVVlU3VRTbhYnn6IpWn7mVwP11IoGgzeCOonKeWgmh4aMQclhHgPGA0kAEeAx4AoACnly0IIAbyAV+lXDMyQUm6o7rxnnNlPtrvxOQ7llYbFTpPHg8dkon3eEVqXFITlnOHEZhGUuSQxVjMuj4cyl8QswGIW3jyvEDhd0n9DT3BYKfFFEQ6rmeJyd0DqzOo7n3bTB/w3Tb0DMJKOqgwjN30j6abgG3plN/hwpeBU2kuhiAB1cFCRVPFdK6VMklJGSSk7Silfl1K+LKV82fe+lFLeKaXsJqXsa8Q5AbSOsbL2wbFhs9Nj8n4FudGtwioiCBcut9eZDO4Sz4d3DOfO87phjzJjMZmwmM1+ZwPeKCulvTdd6bCamTa8iz+dBd7UWZlLkuCwcsfobkxISWJCSpL/xlxY4kTgPceiNRnkFjr9NR3wOoxpb65neXqgjiV4e21TWfpzPbtiG0uCnJFma7DTC1cKrqrPqFAo6p+IpvgiwaBBg+SGDRtYvC6Thz4KuQa4xmhRVIe8I8TXQxQVSlGnRUpV7RtjNeNyuylzn3zdwmomRydCiNGl7DTHszz9EAu/243TDef2OJm+04sKwJtOA69zK3K6kb7nS24bBpyMNmKsZuxWM1N8N/WXVu/GI/Ffsy6pLO0aWqpSf87KBBDhREVUCkWYaYgRVKSZmtqZuZenEGevuxBRi6Jy6imKCr6GWcAtI88IWQMz6TYWOd34fAjgjaoWzRjC6B4ni/3FTjcxPuGCFg0s/G6336kB/ghBLyoAr3OzWQSFPuekOUctmtCiDZfHQ26hkxe/2c2L3+z2K+K0KKSySEdPZZGKdo1RyQkBkU1VAohwooQECkXDodFGUHrOenwFJ0pcdTpvfUdRwQhOChwSY61kFzqJMhHgWGwWwbAz2vDdzhx/dKHVivq2b+kXQ8RYzSy9bZg/GtBHYdo1NBVbsKhAH8G0iKoYjQFMfmWt36lp59RShkbRzqGlIqurWekFFtPCEEHVpzhCoWjWNMcISs/Gx8bTso6RVH1HURqCk5GKW3pFDC7pveF3aRMkp5aSLVn5AdGFVu9Jz8onxmoGvFHU8vRD/migh07irQVkHgnpWflM8601WrQmI+CYO0Z3Y/rwLiGjienDu/ijrRiruVLnZKSeI6m+ZqU5k+nDu7DktmFhcShG6mQ1RdWvFIrw0iQcFMADE+ouizZ5PJRF2Thhjw2DRcaQBKb8ylyS3EInC7/bzc7swoB9y9z4nZEW0ZT6opDJgzsxuEu8/5xL0/b7U205xSdrVG7prSvpHZyWOtNqPJo9laXqJqQksfS2YXx4x3CWVuEwqnICg7vE+8UffX3CjhPFTia/srbCDT4SziQSqbxI2KlQNGeajIOamtq5zuuY/FFUzKlR9OnrTWVuiDJ5I552sdaAfYLrMnarmQkpSaRlHPfvF7zQ02E1YzOfVPdpjmfK4E7+701Q8SZrJCpYnn6Iya+sZfIra3l2xTb//lU5gS1Z+f4obktWPgAuT+hFvJFwJkbqZDVF1a8UivDSJGpQGvcs+YWPN2bV6fzC40GaTHQ4cYT40vDVonrqFsiGQq/iE4DV7HVS+jZC+rVLWt0pPSvff0N88Zvd3nOZ4cPfjzBsm151BwQo8KpStWmpt1Kd0CJYfRe8r1bzCb7mojUZuDweLCZTRFV6CoWinmmIC3UjRVUOSmPEUys5cKJuC3nt5WV0y90fke4SsTav+EATQQTLx7VOCd/uyMFqhmHdEvyOqDKnoe95p5eG15VnV2zju505JLd1kFPsDBAVaNd0WM1+1d+gLvEBtoLXOQVL0RUKRTOhDg6qUTSLrSkhlhPVCOHxUBplI8/uoFVp5VFPTbCYvCksgIIyt1+Rl56VT4nTHaDWS2nfki1Z+d6FtW4Cak7gTSXpow9tm7bWaZCvFhUOtFTczuxCPPJkbSvYjqoiniVp+ytI0RsztVEAKtWgQlFzmkwNSs+sMd3rdLyMQC3KFTTRo8wlWbsnl4ISZ4BkGyAt4zh927dE4E3XBd/UQ9VPJqQk4WjhTQem+2o64aCydUmV2VHZOWKsZlpEmcNm16mkNmIIJaBQKGpOk4ygpqZ25tFP0is4hZrgjaLsYY2iIDCSCtU5ArypsrW7cwLUdBraL/G+7VuSlnE8oG4zZXAnFq3JoMQnM9dqPXX55a61Q6oLE1KSAm7Q+tRfY4wqQkWwkThGoWjuNMkICuCJy1JIirNjqWURSYuish2tCec4wxHdE7jzvG4kOKyM7pGALSioSHBYmTa8i79jRJmbkH3uvtuZQ5HTTZlLUuR089aaDBatyaA4qNtCTX65R3IdTyiFW12iilO55qg2CsBIqAYViqZOk3VQU1M7s/bBsTwxKSVAvl0TTB4PZRYrx6LjwmbX2j25/ptVnw5xOFp4ZyXp03kTUpL8i2vNwqvge2tNBtPeXE/f9i39KTf9YllthpIksNZTE+mz1iD2Ld96qHA6gVA36LrIssOdMlOLbBWKhkeTdVAaU1M7k9iy+pHgodDWRWXHxOMSNf+qQvnFMpdkefohv7Itt9DJgRMlCOGNlrQbrra41iPxO6DcQidpvtZEfTrEBSyWnT68Cw6rOaCzg5ZC0xyA0Ruwllo00uGhLjf1ukQV4V5zpGpECkXDo8k7KPCKJmrbCkl4PLjNFnJiWtX4WO1GH7zYdmnafhZ+t9s/p6nIeXJm04liJ1f+8wcKS8v952jhG3+e4LBSWu4mt9DJy6t3V3AONt2Y9Glvrmfh93v8EZGRG7B2De0c1TmB+rqph3KE4U6ZqUW2CkXDo0mug6qM5Ic/p9zIWNoQCOkhOWcfVnfNm9LazGAxe9c6WcxeJ6Itqq0Ki8kbQekn3F664IcAZaG2rki/NkqLtjRifI6rNuMqQgkZ9EKN4DVP4SR4IXBzWUPVWMUjCkVImnuzWKPUVuYsPB6kMJEdU7v1RU63N0oqc0NJuVf9ECyOgIopQZfnZFNXDYfvQJuZgF/8+gigb/uWmIS3e4UWEWk3Ok1QoaeqVJ0WJWk1MO3mqe8JGKmbqHYdAackujlVdSmVblQovDQrB1XbhrLSZAIpOR4dR4nFVuPje/hEEIB/sesto7oFiDd6JjoqHBfshJanH6LAt6K3zB3YZFaf8tIW1+YWO/0NYfU32eAYsqoboub4tKhsqe+XfX04DO06+t6B4caIc65vR6HSjQqFl2bloKamdmbSWe1rebT3tn7E0drw4l2zgDvP68aBEyVIvP32tBvPhJQk7hh90knlFjvxTcvALLyO6ZZR3QJuzEt0N0pN3VeVU5nsWxeljW4PrjGF2j8YzfHpR2/Ul2S6Pq5jxDnXt6NQknSFwkuTXKhbFc9PGQBQ86aywoSQHgrtMRTaooktK672ELf03gBdbm/UYzGZAmooE1KS2Howj9U7cigoLQ/ouQcnHVKo1kJApQs/9YtrtbqToPJFt9p2LZrQN3TV10Ka4g2zqgW0wT8MwvH5VX1JoTBOs4qgNJ6fMqBWa6OkL1F3JMZ4FKWJFbTZR8EpJW3URJlLBtR09JGPhv6Xtf55ZWkq7bU2YqM6gqOJ5lALqS5aCfd30By+U4UiXBhyUEKIc4QQU4UQN2qPSBsWaSb2r0WqTwhvCySrneMGhhrG2ryDAcFbe0rLOF5hMaw2rM9mESF/xev9aGWOqLKbXvC8qGCCzxec0lK1kPB/B+o7VSiMU63MXAjxDtAN2AhoXU2llHJWhG0LSV1k5sH0m7OC/NKay8YBolzlJOfuw1TF92cSsGzmCCa/spYi3+RbreNDjNXM0tuGVRhZAfhrRMGdwiubzaSfrRTct6+y7c1Rvq1QKE4BER63MQjoLRvbgikDjOmVWKsBhyaPh3JLFDnRcSQWnah0v+S2Dqa9uZ6OrVqw/WghpS43ZnFyXRJ4f1G/tSaDQl1H86Vp+0Omnfq2b8l3O3NI8UVdGsH1IX2dI5Tj0SIu/ej3hkhjrtc0ZtsVioaCkRRfOnBapA05FazadrRWx2ktkHJj4nGZKi5oMgtvam5XdiG5hU52+Cbpujz45z5p0u8JKUnYrCfP4bCaSWnfMmQqTz8mvSqqq3NUJt8OTvmd6v50jble05htVygaCkYcVALwqxBihRBimfaItGENHeHx4DaZORpi8a5bekXpWtMKiddpwckFuvqbl+Yw7jyvG0tuG8aWrHy/QGLyK2uZ/MpalqcfMly/qG6/yoQBRkUSVTmucDq1xlyvacy2KxQNBSM1qHNDbZdSfhsRi6ohnDWoxesyWbBqF3eN6c78Vbs4nFfDMfG+7y45Zx82d7l/c7tYK9mFTto6rBwpcAYcYjPDLaO6VTqJdnn6IRZ+v6fCrKi61ImqalcUvE1vV2U1rMpqYdW9p1AomiGRbHXkc0TbgFjf47dT5ZzCjTaSY2pq51pN4RVIEIIjsW0CZOdHCpxEmajgnOBkik/bPzjiWJK2P8A5aWM4tCGENUE7t6Yc1EdCoaKj4MhKew1UqfbToyIHhUIRLqp1UEKIa4D1wNXANcA6IcRVkTasvpma2rnGx0jf4t18u4MCW0zAe2UnNQ9YdN+yzSICnIPmPBZ+t9s/78lhNWMze8UUvz+vG44W1oAhhEapqpedEUdSmYOrbu1Qk1PTKBSKU4IRFd/DwGAp5VEAIURb4Gvgg0gadiqIMosadzvXFu8ejm2Dw1kcUnauHz1vMZn8XcBT2rfk2x05gLehrNaAdcltwyqcI7jbQWUqMf12fZeEYGei7zShtUGqrCYVUwO1X6jR7k0NpdBTKOoHIyIJk+acfOQaPA4hxIVCiO1CiF1CiAdCvH+6EOIbIcQvQojNQoiLDNodER6f2Icocw1bTPgW7zotVrKjQ3c7104p8K6BSss4zqIZQ/zDBwHO7ZHgdwLBab9QEUuoLuP67ZqDqK5LQlHQiHg9WpQ1vQbNWptDik8p9BSK+sGIo/nCp+CbLoSYDnwGfF7dQUIIM/AiMAHoDVwrhOgdtNsjwL+llAOAKcA/a2J8JGjjsGG31MxJad3OcxytKDVHAV5npLVT0qv5AMpc7oDjHVYz94/v5XcCRm6AobqM67cbcRBTBnfyT+GtqlFsTaKE5tDotDk4YYWiIVBtik9KOVsIcSUw3LfpVSnlRwbOPQTYJaXcAyCEWAJcBvyqPz2grTqNA2q+ajaMaEq+WrTpQyCRwsTh2AQ6n/BGM5UJJLWU3/ThXUI2Kq2qgamGvsGrft+aNHWtbl+VygpNU22cq1A0NCI2UdcnpLhQSnmL7/UNQKqUcqZunyTgSyAeiAHGSSl/CnGuW4FbAU4//fSzMzMzI2KzJjtP7dqaZZuy/GPYjSKkd7BhxxOHaVVaWOl+PRMd5BQ7G/yNX0nGFQpFnYmEzFwI8YPvb4EQIl/3KBBCVN3KwDjXAouklB2Bi4B3hBAVbJJSviqlHCSlHNS2bdswXboimuz8+SkDeHJSSo0jKX+3c0cb3BU/BuBV8R04UeJfiKvVjxpKFwf9dVUqS6FQnEoqTfFJKUf4/lbftjs0BwH9na2jb5uem4ELfddZK4Sw4+1cUbseRGFEk50//FG6cdm0TzBRbokiOyae0wpzK+xS5pKU4W0W63K7yS10s2hNBnarOUDcYEQNtzz9kH8cRygVXmVUlboLkMA38VqSQqFo2BhZB9VNCGHzPR8thJglhGhl4NxpQLIQoqsQwopXBBHcImkfMNZ37jMBO5Bdkw8QSZ5avq3Ga3r8gomYVpRYrJXuV+x0+2tRLre7VqMuqlPhVXVcQ5siq1AoFMEYWQf1ITBICNEdeBX4BFiMNyVXKVJKlxBiJrACMANvSCm3CiGeADZIKZcBfwBeE0Lci1cwMb0hdU13BqntjCPBJ5jocjzL2w3CIkDKgE4SmrovynxyXpN+emt10YvWCV1CjRxKfU+RVSgUitpgpBffz1LKgUKI2UCplHKBEOIXnzS83glnL77q0OZF1WYBryaY6HDiCPGlBd6UnsdToceeScAdo7sB8NLq3Xhk3fruhQMljlAoFGEjkr34gHIhxLXANOBT37aoml6oMfLAhF4kxdl5fGIf5l6eUqNjpU8kcSS2DTab9+vSO6eeiQ4SHFaS2zr45ze7efEbr3MyidDRUH2KJsKR5jvVozoUCkXjx4iDmgEMA/4qpdwrhOgKvBNZsxoG+maytenVJzweXGYLGdY4ipxu9E0qDpwoYdGMIew4WhhQ57pjdLeQabX67F4QjsW2qtuCQqGoK0a6mf8qpZwlpXxPCBEPxEopn64H2xoUi9fVfO2VJpjIjY6jOMqGW55se+TyeJj25vqQx4WKPhqbeKGx2atQKBoe1YokhBCrgYm+fX8Cjgoh1kgp74uwbQ2K+at21fJIr2DiUGwCXY8dBOntUg7eFkXBLE3bH9C+SD/6AhqPeKE23RZU5wqFQqHHSIovTkqZD1wBvC2lTAXGRdashsesMd39vfWM0qqFBYQJ4fFQYm1Bjm/6brHTzeAu8f7punpyCp0kRFtDRh/Vpc0ae91HpQUVCoUeIw7K4mtJdA0nRRLNjqmpnXlyUgpJcXY6trID+P9WxokSF+BL9QHZjtaUWGxIIC3jeMDMKD3bjxb6a0DBnR0cVjMFJU7/GHg9kbzB14fzU2lBhUKhx4iDegLvWqbdUso0IcQZwM7ImtUwmZrambvGdOfACe9o+AMnSgOGEQYT0BRdepBCkNWyLR68Yzf0jO6R4H9u0x0Y3FHCZjVT5ibk4ty+7VtiEpDSviXhpj6im+bQCV2hUBjHiEjifSllPynlHb7Xe6SUV0betIaJvhZljzLh8lBp6i9gyZM/1WcnO+bk3Cj9sXee140Eh5VbRp7h36aPKpanH6LU6cZm9jqx4DHwW7Ly8UhIzwpXq8STqOhGoVDUN0ZaHfUQQqwUQqT7XvcTQjwSedMaJrPGdCcpzs7cy1Po1c7bptBo13NN1ael+vTHfrsjh0VrMigJiqz0UYXW2ii2hRWHPapCFBVJJ6KiG4VCUd8YSfG9BjwIlANIKTfj7avXLNGvjdp8MK/Gx1sEIARZcd5Un0l4H1YzVfbV06Inh2+4oOaMUtq39NeGQjmRxi6cUCgUzRcjDipaShm8YMcVCWMaGxP7t8ckqNEEXhfejuclUXayY1rTIsrMspkjuGVUN2xm7yTeUDUkLXqyW81+CffkwZ34bmdOlbUhpYxTKBSNFSMOKkcI0Q3fxHLfIEL1cxx4fsoA9sy7mEcv7VOj406m+uIpsti4ZMEPfLBhP06390vekHG8wjGh0ndL0vaHbI+kZjopFIqmgJFu5nfi7WLeSwhxENgLXBdRqxohLe0WnC43pa7KC1Imoa9XeRfw7mrRmjOKSzhScHLRbrnbzeRX1gIn5zyFWviq70quf0/NdFIoFE2BKh2Ub7rtICnlOCFEDGCSUhbUj2mNh/mrdpFf6qKl3UKp62T2M85uQYiT66E8IVR9pVF2sh3xtCs8hsNq9vfl02ToVQ0sDOW0tFpVjK9WpVAoFI2VKlN8UkoP8Eff8yLlnEKjKfuC50fllbrIK6m8XOdP9cXEk5AYz5LbhrH0tmFMH94Fh9VcKyej1apaWM0qclIoFI0aIzWor4UQ9wshOgkhWmuPiFvWyCgqc4VM71WvQJcgBGvdDuZ9sc3fj27a8C4svW1YyAipKlVeJGpOSgmoUChOBUYGFu4NsVlKKc8IsT3i1OfAQqMMnbeSw3ne7hImAQ6bhfzSipHTWR3j2HQgr4LTEh4P0mQiofA47Yty/anAGKvZX4PSOBXDBNUAQ4VCUWsiObBQStk1xOOUOKeGyqwx3YmzW2hpt/DkpBT/oMPgNkh7c4oCZkJpaKm+HEc8eVHR/u1FTjcvrd7Nsyu2GVLlRSrSUUpAhUJxKjASQdmB3wMj8GasvgdellKWRt68ijTECKoytJHxdosgPsbG8aKyylV+0gPChMXtosex/UR53GhT5jX1nzYevrLaUqhIp6GOsGiodjUE1HejaFJEeOT720AfYAHwgu95s5ioW1e0SOrRS/tw15juIZ2TFlDZLCZM0juBd1/LRH8fP5OAUckJfif10urdNao/NdSFug3VroaA+m4UCi9G1kGlSCl7615/I4T4NVIGNSX0o+KHzlsZch+Jt/2RMJnwCBDSQ6EthlxHPAmFx4kyC/p0iANg9Y4cPPKk9Dz4l3Z1a6UaEg3VroaA+m4UCi9GUnz/Al6QUv7oe50K3CmlvLEe7KtAY0rx6Vm8LpOHPkqv9P2Wdq+wIsoE5b7cXpfjWTicJSQ4rP4puwKItpoZ3CWe1TtyAK+YYultwyqcU6WKFArFKSfCKb6zgf8JITKEEBnAWmCwEGKLEGJzTS/YnFi8LpOh81ayeF0mU1M7M/fyFKJ8KgmLODnw0CLgjIQYTAIu7tfeOw9KCA62TMRmiwpoDhttNVPkdPudE5xMEwYTnCpScnGFQtGYMOKgLgS6Auf6Hl192y4BLo2caY2f+at2cTivlAW+GVJTUzvTxuEds9G2pZ0fHhjLaXF2XBI2H8zDI2Hd3mNYLWaEx0O5JYroHt0Y3+c0wJsOHNwl3j9DSgAOq5lpw7uEvP6UwZ2IsZr9c6MWrckgt9DJW2syIvvBFQqFIgwYcVDJUspM/QMYrXuuCMHidZkUl7mIs1u4a0x3/7YTxd6ee8eLnSxel8msMd2xWwQe6e2KfteY7jx40Zm0jWtBtEXwfbaLV7YV+qOh9Kx8v2ji3B4JLAmxmFdjQkoSdl/EpS+4GxlfpaIthUJxqjHioB4VQrwkhIgRQrQTQvwXFTlVi9afT/qeL16XyfxVuygt9wBQWu5hwapdTE3tjNNXc9L+zl+1i3vGJXP14NMBeGZLISI6xj/OXT8599kV25j4wg88u2JbSDv0yr7pw7uQ4LAyvZKIS49SkikUilONEZGEAP4A3Obb9KiU8r1IG1YZjUUksXhdJgtW7aKozEV+qYukODt3jenOE5/+Smm5B3uUiQv7nMaPe49xWqyNzQfzmNi/PT/uPcbhvFKS4uxI4PCJYhAmrC4nZ+Qe4LQYC5N1Kq+XVu/2r5FaNnNE2Oxfnn4oZKd0hUKhqBERFknEA0OA3UAZ0NnntBRVoE3e1dZC3TWmO1NTO9Mq2gpAfLTV74yOFJSxZ97FPD9lgL/xbGrX1uTkl/q7njstVnJat6NPUkuW6ByHlu4blZwQ0o7apuoqG/GuUn8KhaK+MBJB7QCeklK+IYRoATyNdwTHOfVhYDCNJYKqDC2y0upSTy33pubOSIhh04E8bBbBo5f28QssNExIPAja52fTujjPcF+8cPfRU335FApFjYhwBDVOSvkGgJSyREo5C3igphdSnET7SbB+7zHyS70pwI2+JrKlLslDH6UztGtgw/gpQ7z1qEOxbSix2klp3zJkJBMc4YS7j1599OVTUZpCoQBjDmq/EOJ6IcSjAEKI0wFDffiEEBcKIbYLIXYJIUI6NSHENUKIX4UQW4UQi42b3jjRS8+XbcqqdL91e48x9/IUkuLszL08hblX9GP4Ga2RwsTRhPakHSoKKWIIFjdUlqqrLeE+XyiUQEOhUIAxB/VPYBhwre91AfBidQcJIcy+/SYAvYFrhRC9g/ZJBh4Ehksp+wD3GDe9caLVmO4a052J/dtjEt4xHElxdiad1d7fFV1LARaVuXhq+TYWr8vkrZtTSU5oQYHHxB5HIvHpVgs/AAAgAElEQVQOW4VIpil0Hm8Kn0GhUNQdIzWon6WUA4UQv0gpB/i2bZJS9q/muGHAHCnleN/rBwGklPN0+zwD7JBSLjRqcGOvQQWjyc9n+UQUevRzpgQQa7cgpaSwpBxpMnFxeysvnBOP0qwoFIoGS4RrUOW+aEgCCCHaAh4Dx3UA9DmaA75tenoAPYQQa4QQPwohLjRw3ibFU8u3cTivlKd9Ygl9eyRtzpTA++Xnl7ood3vAZCJKwGdZTl74tfCU2q9QKBSRwoiDmg98BCQKIf4K/ADMDdP1LUAyMBpvCvE1IUSr4J2EELcKITYIITZkZ2eH6dINCy2O1deopqZ2ZtOc8fz18hTsFuHvuScBR4soBPD3X4v44mDtRnMpMYJCoWjIGJmo+y7wR2AecAiYJKV838C5DwL6IkJH3zY9B4BlUspyKeVeYAdehxVsw6tSykFSykFt27Y1cOnGg7ZO6oEJvfztkVraLaR2bR3QaLZVjA0J2CxmkuLszB7fkz+O7wHAnT/ksnBD5U6mMkekxAgKhaIhYySCQkq5TUr5opTyBSnlbwbPnQYkCyG6CiGswBRgWdA+H+ONnhBCJOBN+e0xeP4mgbagd2pqZ397pBibxb+IV2s0q4kr/jShl3//20d3xyrAbTLz7M5yjpWFzrxW5ojqS4ygIjWFQlEbDDmo2iCldAEzgRXAb8C/pZRbhRBPCCEm+nZbAeT6BiB+A8yWUuZGyqaGjl7hp38OgY5MQwiBzertfF5qsTLly0OUeyqKXipzRPUhGQcVqSkUitphZKJurZFSfg58HrTtUd1zCdzneyg4WYvST+PVE6z6e/CiM3n4o3SE9LCj1MwTP+fxl0GBZbxQk3brEzUhVqFQ1AZDDkoI0Q4Y7Hu5Xkp5NHImNV+CBRKL12X6WyE9MKGXPw14OK+URz72TufVnNjfv9zBiWIn7+wtxVWcRdbuAw1mku6pdpAKhaJxUm2KTwhxDbAeuBq4BlgnhLgq0oY1R4LTelpNKr/UFVCLMgnwSAIGIf705/N5+sp+ACw5LMhwmlm0JkPVfhQKRaPFSAT1MDBYi5p866C+Bj6IpGHNkeC03qwx3Xl6+TYkBNSiwLt+6nixk16PfI7VYvZHWDuPFPDK93s50Oo0HAWHKSosYWnafhXBKBSKRocRkYQpKKWXa/A4RR3R1kFtnjM+wHFNTe1MtM1CabmHUpckv9TF08u3MXTeSjq1ieby/km4TWYy45OIjo02VPupidJOqfIUCkV9YMTRfCGEWCGEmC6EmA58RpDwQVH/aF0m7BZBnN3iHW6YV8qL3+zmmWvO4tzubSj0mMhK6MDZ3RKrPV9NlHZKladQKOoDIwt1ZwOvAP18j1ellH+KtGGK6mlhs3BhShItbBbG9Er016/e37Cf344W0tpu5kCJ5Iovj5BfXnV3qpqsiVLNXBUKRX1QZbNYXw++r6WU59WfSVXT1JrFGkFT8zldbn+9SVPzaYKJpDg7ax8cC5xsMisAPB6kyURqGzNvnZuA3awayyoUinokUs1ipZRuwCOEiKu1cYpaoW8aq6n5tHrTIx97BxomxdmZ2L99gPIPTqoB+3eMA5OJFmbBulw3I5fu4bMtxutG1dWaVC1KoVBEEiM1qEJgixDidSHEfO0RacOaO/o1Ufp6k8AbMX2z7ShFZS5WbTvKXUGjOrSuE4cLypBATIsoTEiyo2L46+YCqhuxolFdrUnVohQKRSQx4qD+A/wZ+A74SfdQRBD9mihNzbftyYv4q2/KrjZ+Q79GCiqO60iKs3Pf+T24eeQZCCnJsrbkH1uNjeiortakalEKhSKSVDuwEEAI0QI4XUq5PfImVU1zrEGFYvG6TP8aKW0NFJysP+lrUtr+z6zYTl5xORL468BYrusWc2qMVygUzYc61KCqXagrhLgUeBawAl2FEGcBT0gpJ1Z9pCKSVNarb9aY7ixYtSugJgXelOGJ4nLi7BbySl38+ecCHFEmLju9RX2ZrFAoFDXCSIpvDjAEOAEgpdwInBFBmxR1YGpqZ+4a0535q3axeF2mf7t+XMf95yfjAe5dl8cn+0r8+yjRg0KhaEgYGvkupcwL2mZk5LsiQujrTKG26QUWGvpxHTPH9uDuMd0rOCklelAoFA0JIw5qqxBiKmAWQiQLIRYA/4uwXYoqCOWAglV/wdJzDc2RtYuzBzqpzGIlelAoFA0KIw7qLqAPUAYsBvKAeyJplKJqQjkgTYpeVOYCCJnmg0BHdu8FPU86qfX5OGPj6mWAoUKhUBihWhWfEGIk8D/fol1t20Ap5c+RNi4USsVXOf3mrCC/1EWc3UILm6VSNV/wjKnnvtzOP1btwgQ8N6Qll3WOPkWfQKFQNDki1UnCxwpglRBC33F0YU0vpIgsi9dlkl/qjZ4koaMsONkJXb9+6t4LenLP2JOR1MeZxdVeTwkqFApFpDHioLYDfwO+FUKc49umGro1MOb7nI1JnIyMNFFEMKGcV0ZuMQIPHuA+A05KCSoUCkWkMeKgpJTyU2Ai8IIQYibeH+mKBoTmdJ6clBLSKQUT/B9w2aYsJCbDTkoJKhQKRaQx4qAEgJRyJzASGIV37IaiAVFVxBSMJpR45ON0v4hiYv/2mARcdlZHXbovjwveTven8fRpvQkpSRUEFc017ddcP7dCEWmMzIMaoHteJKW8BrVQt1Eza0x3/5gOrQ71/JQB7Jl3Mc9PGcA95/fkvnHJSAQ7WiTwt40nkFLWuHlsc7lxq3SnQhEZajW6XUq5L9yGKCJH8MLeqamdeXJSSqVrpQBmjevBJX1PAynZY4vn8V/ySUlqiUlASvuWIY8JTvs1lxu3SncqFJHBULPYhoSSmYdG6yAxK2j0BpxsIKvJz0PtUxmfbs7iviUbcXokic5CEo4dJtFhZdGMIdUeuzz9EEvT9jN5cCe1tkqhaK5EQmYuhLjb93d4XWxT1A+huktoaAIKCRVqT9VxSb/2LLppCA6riaNWB4fbtGfiwOojheXph1iinJNCoagDVaX4Zvj+LqgPQxR1o6r2RpqA4oEJvSrUnoxwTvcEltx2DgnRURyLiub1w2aOlrqrPKa5pPcUCkXkqMpB/SaE2An0FEJs1j22CCE215eBCmMYUfEZqT1VRkqHOP5z5wi6xNvZmu/hqq9zyCh0+d8PFkSouoxCoagrVdaghBCn4e0kUWH2k5TSWI4ozKgaVN0J1e7IKDmFZcx4/Ue2HCok1uQhufAIvzv7NH/ElGCwPqVQKJoJkWp1JKU8LKXsDxwCYn2PrFPlnBThYf6qXf5x8U8v31ZhdEdVJDhsvHf7cEae0ZoCj4lN9kT+sf6IipgUCkXYqVZmLoQ4F9gJvAj8E9ghhBgVacMUkUPrfN7SbvELJ4w4Kk2uvmzjQV6/KRWrALfJzDZHOw5GxfLm9MEVBBHNZS2UQqEIP0bWQf0fcIGU8lwp5ShgPPCckZMLIS4UQmwXQuwSQjxQxX5XCiGkEKLGIaCi5kxN7cymOePZPGc8D0zoFaDwq0o8oVcKWi0mJvRLQiBBCJ7cXMjUFYcocwemjJVYQqFQ1BYjDipKSrldeyGl3AFEVXeQEMKMN+qaAPQGrhVC9A6xXyxwN7DOqNGK8KFX+FUnntCUgqldWzN03kq+2Z6NRGi9sFhbYGLKNzkcLTmp8AuV+lNRlUKhMIKReVBv4B3x/i/fpusAs5TypmqOGwbMkVKO971+EEBKOS9ov+eBr4DZwP1SyioVEEokcerRFv62tFuIsVlI7dqaH3bl4Cp3c6LMzWk2eGVEG/q3Dv07Ztqb65WgQqFoLkR4HtQdwK/ALN/jV9+26ugA6PM6B3zb/AghBgKdpJSfGbJWUWuC2x3VBS2SemBCL9Y+OJbnpwxgwyPn8/Xs8xh8ehyHy+Cab3L5ZF9JyOMToq0AtPH9VSgUilAYaRZbJqX8PynlFb7Hc1LKsrpeWAhhwlvf+oOBfW8VQmwQQmzIzs6u66WbHYvXZfLIx+nV1phCHad3atprIOSaqwSHjXdvPYdrB3WkzAN3r8tj3qY83EFR+s7swoC/CoVCEYpaNYs1yEFArznu6NumEQukAKuFEBnAUGBZKKGElPJVKeUgKeWgtm3bRtDkpsn8VbvwSO8ww5os0A1un1RVOyUNq8XE3Cv78cTE3pgFvLKjhFu+P0Z+uce/z6jkBEzC+1ehUCgqI5IOKg1IFkJ0FUJYgSnAMu1NKWWelDJBStlFStkF+BGYWF0NSlFzajrMMPg4zalV1U5JjxCCG8/pyjs3p9LKbuabI+VM/DKbV9IOMfmVtaRlHOeO0d24f3yvOn0uhULRtDEikugrpdxSq5MLcRHwPGAG3pBS/lUI8QSwQUq5LGjf1SiRRJPjhVU7mf/VdpxSYJIe2uXn0Lokn7Y+Zd+StP1MUQ1lFYqmSx1EEhYD+/xTCGEDFgHvSinzjJ5cSvk58HnQtkcr2Xe00fMqGg//WrcPpxTEWEwUueBQXCJl9mhm9I0JWCOlHJRCoQjGiEhiJF5peSfgJyHEYiHE+RG3TNEk0NKCD1/am39eN5BYq4ljNgfP7TMxuI9qj6RQKCrH8MBC38LbScB8IB8QwENSyv9EzryKqBRf42b/sWJm/msDm7IKMAP3p8RwWy8HJiFOtWkKhSISRHIdlBCinxDiOeA3YAxwqZTyTN9zQy2PFAqNTq2jef/3I/jdiC64gafTi5j+bS451cyXUigUzQ8jKr4FwM9AfynlnVLKnwGklFnAI5E0TtE0sVpMPHxJH96YPoh4u4Xvsl1MWJHNmiN1Xl6nUCiaEEYc1MXAYillCXgX2AohogGklO9E0jhF02ZMr3Z8fu8ohnRuRbYTrvvuOI/+dIIil6f6gxUKRZPHiIP6Gmihex3t26ZQ1JmkuBYsvnUY947tjsUEb+8p5cIvsvkx23mqTVMoFKcYIw7KLqX096TxPY+OnEmKpoSRHoAWs4m7z+/JJzNHcGa7GPaXSKasPsacn/MoVtGUQtFsMeKginxNXQEQQpwNhO4Cqmi2VOaIjLRH0ujTPo5P7hrF3WO80dSi3SVMWJHN+kqiKTW2Q6Fo2hhxUPcA7wshvhdC/AAsBWZG1ixFY6MyR2S0PZKG1WLi3gt68vGdI+iVGENmsWTy6mM88UseJS41DFGhaE4YWgclhIgCevpebpdSlkfUqipQ66AaJovXZbJg1S7uGtO9Rv3+qsLp8rBg5Q7+uXo3bgldowVXtvWQttXbHglgadp+JjfwVknL0w+plk6K5ksd1kEZdVDnAF3QtUaSUr5d04uFA+WgGheL12Uyf9UuZtXBcW0+cIL7l/7CjuxiAFoXnaCXLGTxjBr/ez8lqAGNimZNhBfqvgM8C4wABvsejePOoDjlGKlBVSek6NexFf+9exQzR3fDLOBYTCs2xbbnvd1FeAx2QjmVhBp7r1AoqsdIN/PfgN7SaE+kCKMiqMaFkdSfNkI+Kc7O2gfHVnm+XUcLeeyjzazZexyAs1qZeXJQK1LiQ4+XVygUp5gIj3xPB06ruVUKBUxN7Rxy+q6emggpuic6+Netw1hw7QDaOaLYeMLNxK9zefSnE+Q5aydJV2pAhaJhYsRBJQC/CiFWCCGWaY9IG6ZoPhhxYhqL12Uy7KlVFJSWs3L2GG4Z3gUhvAt8xy4/yocZJdQ02FdqQIWiYWLEQc3B28V8LvB33UOhqHf0NS2HzcIjl/bhs7tHMuT0OHKc8Ie0PK5ZlcOmY1ULTUvLyhhy4430v/ZaVn80jwNbv/TXiPYePEjqtGl0nzSJyQ8+iLPce64yp5PJDz5I90mTSJ02jYysLP/55r35Jt0nTaLnFVewYu3ayH0B9UhGVhYp11xTq2MX/fe/ZGVn19v1FE0TI/OgvgUygCjf8zS8zWMVCj+a0OGeJb9U2zmiLoRKB/Y6rSVL7xjO36/uT0K0hbRjbi5bmctd/zvG/iJXyPPYrFZWvfwym957j13/eR9bSSbxMgeAPy1YwL1Tp7Lr44+Jj43l9U8+AeD1Tz4hPjaWXR9/zL1Tp/KnBQsA+HXPHpZ8+SVb//1vvliwgN8/9RRud/Puzl4bB6VQBFPtRF0hxO+AW4HWQDegA/AyUHU1W9Gs0CKbZZuy8EhYsGpX2NZD6Zma2jnkeYUQXHl2R8b1bsdL3+zkjTUZ/Pegky+ycrixWwtm9o4l3mYK2N8R7e3YVe5yUe5yIYRASsmqtDQWP/kkANMuuYQ5r77KHVddxSfffsucW28F4KqxY5n5zDNIKfnk22+ZcsEF2KxWunboQPdOnVi/dSvD+vULsPH1jz/m6bffppXDQf8ePbBFRfHCn/5E9vHj3D53LvuOHAHg+fvuY/hZZzHnlVfYd/gwew4eZN+RI9xz7bXMmjIFgH99/jnzlyzB6XKR2qcP/3zgAQBu/stf2PDrrwghuGniRO697roAG97/+msef/VVzGYzcQ4H3732Gm63mwdeeIHVP/1EmdPJnVdfzW1XXhlwXFX7PL1oEf9avhyTycSEc85h0JlnsuG337jukUdoYbez9o03+HXvXu577jkKi4tJaNWKRXPmkJSQwE+//cZNTzwBwAVDh9biX4SiKWNk5PudwBBgHYCUcqcQIjGiVikaHbPGdGfBql2kdm3Nur3HDHeOCDdxLaJ44KLe3HBOV/7+xW98tPEQr+8q4d8ZJcw808G05BjsZu9wRLfbzdk33MCu/fu58+qrSU1JIefECVrFxmKxeP/X2Hkc1u/IZHn6IQ4ePUqndu0AsFgsxDkc5OblcfDoUYb27eu3oWNiIgePHg2wKys7m7+8/jo//+tfxMbEMOb22+mfnAzA3c8+y73XXceIs85i3+HDjJ85k98++ACAbZmZfPPyyxQUF9Pzyiu546qr2LV/P0u/+oo1b7xBlMXC7596ineXL6dPt24cPHqU9H//G4ATBQUVvp8nXnuNFS+8QIfERP/7r3/yCXExMaS9/TZlTifDb76ZC4YOReiGSFa2z7aMDD757jvWvfUW0XY7x/LyaB0Xxwv//jfP3nMPg3r3ptzl4q6//Y1P/v532sbHs/TLL3n4xRd547HHmPH447zwxz8yauBAZv/jH2H5N6BoOhhxUGVSSqf2j1UIYQEahORc0XCoLLI5VXRo1YL/mzKQm0fl8dSnW/l+z3HmbSnk7Z1F/KFfSyadbsdsNrNx8WJOFBRw+f33k75rF6clJAScZ9mmLFweWWcBxfqtWzl34EBax8UBcPW4cezI9KZBv16/nl/37vXvm19URGGxd1HyxcOHY7NasVmtJMbHcyQ3l5Xr1/PTb78x+MYbASgpLSUxPp5LR41iz8GD3PXMM1w8YkTIiGR4//5MnzOHa84/nyvOOw+AL3/8kc27dvHBqlUA5BUWsnP/fnqcfrr/uMr2+Xr9emZceinRdjuA//Pp2Z6RQfru3Zx/552A94dBUkICJwoKOFFQwKiB3lafN1x0EcvXrKntV6xoghhxUN8KIR4CWgghzgd+D/w3smYpFOGhT/s43rn1HL7bkc3cT7ey7WgR963PY+G2Au7r25KxSTZaxcZy3qBBfLF2LX+4/npOFBTgcrmwWCycc7qNtbHxTB7ciUM/JbL/yBE6tmuHy+Uir7CQNnFxdEj0btc4cPQoHRKNJxk8Hg8/vvkmdputwns2q9X/3Gwy4XK7kVIy7ZJLmDezYkvMTe+9x4q1a3n5ww/591df8cZjjwW8//JDD7EuPZ3PfviBs2+4gZ/eeQcpJQtmz2b8sGEB++pFIJXtY0QQIoE+Z5zB2jffDNgeKsJTKPQYUfE9AGQDW4DbgM9Rk3QVjYxRPdry2T3n8uzV/UmKtbLl8HFuWnmAS7/M5r978vhq3ToKRUumL0rjzOQUPli5EoAdv67l4akTmZCSxMRRo3jr008B+GDlSsYMHowQgomjRrHkyy8pczrZe/AgO/fvZ0ifPgHXH9y7N9/+/DPH8/NxuVx86ItEwFt7WbB0qf/1xu3bq/wsY4cM4YOVKzl67BgAx/LyyDx0iJwTJ/B4PFw5dixP3nEHP4c4z+4DB0hNSeGJ22+nbXw8+48cYfywYbz0wQeUu7yCkh2ZmRSVBA4sqGyf81NTefO//6W4tNRvC0BsTAwFviiwZ+fOZB8/ztrNmwFvzW/r7t20io2lVWwsP2zcCMC7y5dX+bkVzY9qIygppQd4zfdQKBotZpPgqrM7ckm/JJ56dwXPPDiLL91uvpQeOvcbxXrPGYjCPBJ7jef/3n2XR156iQE9e3LzZZcBcPNll3HDo4/SfdIkWrdsyZK5cwHo060b14wbR++rr8ZiNvPiH/+I2WwOuHaHxEQemjGDIdOm0bplS3p16UKcwwHA/NmzufPpp+k3ZQout5tRAwbw8kMPVfo5ep9xBk/ecQcXzJyJx+MhymLhxT/9iRY2GzMef9zf/mmeL6WmZ/Y//sHOffuQUjJ2yBD69+hBv+RkMg4dYuB11yGlpG18PB//PXAlyS2TJoXc58JzzmHjjh0MuuEGrFFRXDR8OHPvvJPpl1zC7XPn+kUSHzz9NLOefZa8wkJcbjf3XHstfbp1483HHuOmJ55ACMEFqam1/4+raJIYaXW0lxA1JynlGZEyqipUqyNFuCgtd/Peun28vHonRwq9a51iXaVM6WzjwWGnYdKJBIxSVefywuJiHNHRuFwuLp89m5smTuRyXx1IoWiy1KHVkZEalP6kduBqvJJzhaJRY48yM2NEV65NPZ2l6/fx0je7OFwIrx2E777I5q6UWCZ0tGOugaPSd6UIdlBzXn2Vr9evp7SsjAuGDmXS6NFh/kQKRdPC0LiNCgcJ8ZOU8uwI2FMtKoJShIvgUSBlLjf/TtvPS6t2klXgneIb7XZyaXsLj51zGtGW6ku2y9MP1XhGlZoXpWjSRHIelH7cO15RxSDgDill/5peLBwoB6UIF5V1US9zuflgwwEe+3gLLrzRUysLXNctmmnJMSS2MFd2ylqh5kUpmjQR7mau7783DzgbUA2zFI2eyrqo2yxmrhvamTmTUmgVHcXpcTZOuODF7cUM/yybP6w7zra88A2VDjUv6pdt27jZ12EhXFTVS1DPTY8/TuL551foi7dpxw6GzZhB38mTufTee8kvLARgy65dTJ8zp0a2fPG//9HziivoPmkSTy1aVCN7v/rxR86+/nr6Tp7M2ddfz6q0tBpdW9F4MNKL7zzd43wp5e+klFXrYBWKRkB1XdSvH9qFjY9ewHcPjuPDO4Zx4ZltcUn4cF8ZF36Zyw3f5PDt4bIad08PZkJKEotmDAlI7819801/W6OakpGVxWhfSyY9lfUSDGb6pZfyRYj3bnnySZ6aOZMtS5dy+ejR/O2ddwDo2707B44eZd/hwxXPNWcOq4MyHm63mzuffprl8+fz6/vv896KFfy6Z49hexNateK/zz3HlqVLeWvOHG549NHqvxRFo8TIRN37qnrUh5EKRbh59NFHef755/2vH374Yf5RRaudszu35uVpQ1g9ezTThp5OC4uJ73NcTPv+OGOXH2Xh9kJ2HcnlytmzGXzjjQy+8UbW+Nb33P3sszzxmneVxoq1axn1u9/h8XiYPmcOt8+dy6AbbqDHFVfw6fffA1BQVMTmnTvp36NHBTvyCgvpecUVbM/IAODahx7itY8+MvSZP/n2W6Zdcgng7SW4cv36kM511MCBtG7ZssL2HZmZ/q4P56emBqzlunTkSJasWGHIjvVbt9K9UyfO6NgRa1QUUy64gE++/dawvQN69aJ927aAV+JfUlZGmdNp6NqKxoWRFN8g4A68TWI7ALcDA4FY30OhaHTcdNNNvP3224C3k8OSJUu4/vrrK+w3cuRIzjrrLP/jsrHD+WTODcwbCrPH96CdI4o9RZInNxcycPZcLGdN5OXnFvLB009zi6/h7LyZM1n61Vd8s2EDs/72N9587DFMJu//ehmHDrH+rbf47PnnuX3ePErLytjw22+kdOsW0u44h4MX/vhHpj/+OEtWrOB4QQG/u/xyQ5+5sl6CRunTrZvfkbz/9dcB3TMGnXkm3/scck3sgNC9C43a++HKlQzs1Sug44ai6WBEZt4RGCilLAAQQswBPpNSVvy/OQghxIXAPwAzsFBK+VTQ+/cBtwAuvN0qbpJSRmZOg6JREKysixRdunShTZs2/PLLLxw5coQBAwbQpk2bCvt974tqKuO2Ud1Yue0o7/5vL+9mbOSjd/fx0btgNwksZYUcyiskKc7Baw8/zKhbb+W5e++lW8eO/uOvGTcOk8lE8umnc0aHDmzLyOBQTg5t4+Mrveb5Q4fy/sqV3PnMM2xavNi//fL772dvVhbO8nL2HT7MWVOnAnD3lCnMmDixpl9RBd549FFm/e1v/GXhQiaOGoU1Ksr/XmLr1v7xGivWrvWn4/YdPswPGzfiiI7GFhXFurfeqrMdGlt37+ZPCxbw5Ysvhu2cioaFEQfVDtDHz07ftioRQpiBF4HzgQNAmhBimZTyV91uvwCDpJTFQog7gGeAyUaNVzQ99AMJq3JQ4XBkt9xyC4sWLeLw4cPcdNNNIfcZOXIkBSF6xj377LOMGzcOi9nE+D6nMb7PaXz6gJl7X/mQjzcdJbfE2xJo7MpCLuvkwvLbdtrExZGVkxNwHhG0xkoIQQubjVJdymr8zJkcOXaMQWeeycI//xmPx8Nve/cSbbdzvKCAjr4o46NnnwW8Najpc+aw+tVXA86t9QwM7iVolF5duvidwY7MTD774Qf/e6VOJy18vQTHDxvm79k3fc4cpl9yCaMHnRRwGe1dWJW9B44c4fLZs3n78ccDHL6iaWEkxfc2sF4IMccXPa0DjPwMGgLsklLukVI6gSXAZfodpJTfSCmLfS9/xButKZoxlSnrgtE7stpy+eWX88UXX5CWlsb48eND7vP999+zcePGCo9x48ZV2PfC8eOx75wJ5mkAABeBSURBVPya/z08jvnXDqCnKZtiN7yzaR9PvfMvkm96nndWfc+nGzb5j3n/66/xeDzsPnCAPQcP0rNzZ87s2pVd+092T1/xwgtsXLyYhX/+MwDPLV7MmV27svjJJ5nx+OP+/njVUVkvQaNovf88Hg9Pvv46t+tmRu3IzKw0LRnM4N692bl/P3sPHsRZXs6SL79k4qhRhu09UVDAxffcw1MzZzL8rLMM269ofBhR8f0VmAEc9z1mSCnnGjh3B0A/o+CAb1tl3AyE7BYphLhVCLFBCLEhW03pbNJUp6zTMOrIqsJqtXLeeedxzTXXVOidVxvmz5/Phg0bGDxwAA9cO46ux37kq3tHYlv7Cknn38IBc2vc593FFY88weQvs8gsdJHUrh1Dpk1jwqxZvPzgg9htNnp16UJeYSEFRUUVrrE9I4OFH3/M3++5h5EDBjBqwACefP11Q/bdfNll5Obl0X3SJP7v3Xd5ytcNPSs7m4tmzfLvd+1DDzFsxgy2Z2bS8aKLeP3jjwF4b8UKelxxBb2uuor2bdsGpA2/2bCBi0eMMGSHxWLhhdmzGX/XXZx51VVcM24cfXzO7dGXX2aZr85Vmb0vLF3Krv37eWLhQs6aOpWzpk71O09F08JQJwkhxAggWUr5phCiLeCQUu6t5pirgAullLf4Xt8ApEopK8wIEEJcD8wEzpVSllV1XrVQVxEuPB4PAwcO5P333yfZNzwwUqTO/Zoj+WVEm8GFwOmW5Hz2HHHJg5k8dgxXdI1hWKLV31bpuXffJTYmhlsmTYqoXcHUpqtFmdPJubfeyg8LF/oHPSoUfiK5UFcI8RjwJ+BB36Yo4F8Gzn0Q6KR73dG3Lfj844CHgYnVOSeFIlz8+uuvdO/enbFjx0bcOQHcPTaZpDg7j0xMIe2R85l3RV8SYqJweuA/+51c/91xhi47wqM/nWBdtpNbr7wSm06EUF/oewkaZd/hwzw1c6ZyToqwY6TV0UZgAPCzlHKAb9tmKWW/ao6zADuAsXgdUxowVUq5VbfPAOADvJHWTiMGqwhK0ZTIzC3io58P8J+f9rPvxMnfZ22tMKGjnYtOj2ZwQlSNGtbWhdr0ElQoqiTCvfjWSymHCCF+llIOFELEAGurc1C+Yy8CnscrM39DSvlXIcQTwAYp5TIhxNdAX+CQ75B9Usoq9bDKQSmaIlJK0g/m89nmLD7bdJD9eYHO6sIOdi7uXL/OSqEICxF2UPcDyXjl4vOAm4DFUsrQfVIijHJQiqaO3ll9vvlgQGTVJgrOS7IxrkMLRp5mJcZAh3WF4pQSSQcFIIQ4H7gAEMAKKeVXNbcyPCgHpWhOSCnZmpXPp5sqOiurgKEJFsZ1bMHY9nY6RIe3y7pCERYi5aB8i22/llI2mLGfykEpmitSSnYcKeTr3w6zMv0wvxzMDxh13SvWxLgOdsa0t9O/tUoFKhoIkZqoK6V0CyE8Qog4KaXxpl0KhSLsCCHoeVosPU+L5c7zkskpLOObbUdZufUQ3+/KZVuBh23binlhWzFxFhjeNooR7Vswsp2VTjFKYadofBj5V1sIbBFCfAX4Vw5KKWdVfohCoYg0CQ4bVw/qxNWDOlHmcvPjnmOs3HqIb7cfJfNEGZ8fKufzQ965VV2jBSPa2RiZZGdYopXYqKZZu1LTiZsWRhzUf3wPhULRQLFZzJzboy3n9vCOocjMLeL7nTn8sP0Ia/YcY2+xm717S3lnbylmAQNamUlNtDG0nY2z20QZGmffGNCv45qQkqQcViOnUgclhDhdSrlPShm+9sMKhaJe6Nwmhs5tYrh+aGdcbg+bDuTxw85svt92hF8O5rPhuJsNx4t5cXsxFgF940ykJtpJTbQxKCGq0UZYUwZ38q/jgooOKxwop1d/VCqS0NY9+Z5/KKW8MuSO9YwSSSgUdSO/tJz1e46xbk8O63bnkH64EI/uNmACUuJMpCbaODvBxsA2USS2aJwKwUgsPJ725npyC50kOKwsmjEkLOds0kRIJKGXAJ1Rc6sUCkVDpKU9inG92zGut3dMR0FpORsyj7Nudy7rdmez5VABm/M8bM4r4bWdJQB0tAsGtoliYFsbA9tYObOVhShTw1cJTkhJCnuUExylKSJHVQ5KVvJcoVA0IWLtUZzXM5HzeiYCZ1JU5uLnfcdZvyeXnzNy2XggnwOlHg4cdLLsoHdOld0E/eLMDGhr46w2VlLiLXSMNtdofEdjJRJOTxGaqhxUfyFEPt5IqoXvOb7XUkrZMuLWKRSKeifGZmFkcltGJnsFF26PZMeRAn7ed5yf9ubyS+Zx9h4vZf1xN+uPFwPekW7xUdC3lYX/b+/Mo+Oq7jv++c4iyVqsxcZGwXjDJjEJpJCcQDZKIU1DNnpOaEJDIKFpk5aWNM3SJm0OIelpe7rRNmnSJBCWbCTNclKHQ0JzIC60lB2CCWZxbRzLGNmy9tFopJF+/eNe2WMxwiMbad5Iv8857+i+++5793vvjO5v7vJ+99SOOk5dVsdp7Vk6l6QWhdFy5oYZDZSZ1eags+M4LyjplNjUuZRNnUu5OO7T1Zsb46Ff9vHQrj4e2d3Lo88M0Zsvcsf+cEwZrWVZOLU9wyntdWxqz7KpNcu6lrS/ROxUhL+95zjOrOloquO8TSs5b1OYxzIz9vTn2do1wNaufrbu7uORZwY5MDrBln1Ftuw7tOtvfQpe3JJiU1uWTe11vKQty6bWDK11tbly0Jk73EA5jnPMSGJVeyOr2hs5/9QwP2NmdPXleaRrgG17B9jW1c/j3UPsGRyLizAKsOuQb8HOetiwNMPJrVk2tmbZuDTDhqVuuBYzbqAcx5kTJHFiRyMndjTyltMOLSoYGBln27ODbNs7yOPPDLDtmQGe2Jdjb8HYu7/InfuLQP5g+pX1sLElzYbWOja0ZlnfkmZ9S4aVDT6/tdBxA+U4zrzS2pjlrPXLOGv9soNxE5PGL3tHeKp7iKf2DbN97wBPdg+xvWeE7oLRXZjgv3vylBquxjSsbUwFg9Vax7qWDOta0qxr9l7XQsENlOM4VSedEuuWN7FueRNvfOmh+IlJY09fnie7h3hy3xA7uofYuX+YnQdG6M0XeWxokseGJuGZ8cOe15qBNU0pVjelWd2SZU1LhhOb0qxpznD8kpQv0qgR3EA5jpNY0imxelkjq5c1HnyxeIr+kTF29uTYsT/Hzp5hdnYPsaMnx9O9eQaKk3GeaxI43HjVCVYtESc0plnVnOaE5mwIN6U5oTHNSjdgicENlOM4NUlbYx2nr67j9NXth8WbGfuHC+zuHWHXgRF2Hcixu2c4hHvz9IyMs2PE2DFShJ4iUDjs/oygs0GcsCTFixrTHN+UobMxQ2djiuOXpHlRY5r2Ovn81zzgBspxnAWFJFa0NLCipYFXrOl4zvWRsSK7e/Ps6R9hT1+ert4RunpzIdw/Ss/IOLvzxu78BPROAGPPeUZdCjrrReeSFCsb06xYkmZlY4YVDakQbkixYkmKpgXiJb5auIFyHGdR0ViXObjxYzlGxyfY059nT1+eZwdG2TswyrP9Ofb259k7MMrewQKDhQl25Y1dB41YeZrTsKJeHNeQYnlDiuMa0ixvzITz+hTLG9Isj+H6tPfIpuMGynEcp4SGbJqTjmvmpOOaZ0yTKxSD4RoYpXtwlH1DBboH8uwbyLNvcJTuoQLdw2MMTxjDI8aOkQlggunzYaW0pKGjTnTUi2X1KTrq03Q0pFnWkKajPnXwaK8TbfUpWjILf5jRDZTjOM4saarPsGFFMxtWzGzEzIzBfJHuoVH2DxXoGS6Ev0MF9g/m6YlxPbkxDoyMMzQBQ3ljV96ASaA447MhzJW1ZaEtK9rrUrTVifb6NO0NaZ7ty/HY7j7OPqmdc07qoDWborVOtNalaMmqZhaBuIFyHMeZAyTR2piltTHLySvLDydOMTlpDI6OcyA3Rm9ujAPD4W9vrsCB4QK9QwV6cwV6c2P058fpyxcZGZ+kZwx6xgxyU8OMU0ZN0NDBU3vgq3v6DtcFNGdgaQZasimWZhWPFEvrUrTUpcPfrGjJpmjOiuasaMmI5njelJkfI+cGynEcp8qkUqKtsY62xjpOOq6yewrFCfpHxukbGaMvN07/yBh9I+P058e4e3sPDzzdy/rljbQ0ZBnIjzMwWmRgtMhQYYKhIgwVgdHJo9bclA6Grimtg0arKSOasmFxSFNWNGXTNBeOPg83UI7jODVIfSbNyqVpVi5teM61y8/ZMON9E5PG0Og4Q6NFBvLh72A8H5w6zxcYHBknVwjnQ6NFhgvxGJtgeGyS3MRUx82Yqy0D3UA5juMsItIlvbWj3RN4ctLIjQXDNTJWZLgwwUg0YLmxIrnCBLlCkVyM+/RR5uMGynEcx5kVqZRoacjS0pCtKP3RGih/i8xxHMdJJG6gHMdZcHzrnl2c9Te38a17dlVbinMMuIFyHGfB8bnbt/PswCifv317taU4x8CcGihJb5L0hKTtkj5R5nq9pO/E6/dIWjuXehzHWRx86NwNdLY2cMW5M69mc5LPnC2SkJQGvgD8OtAF3Cdps5k9VpLs/UCfmW2QdBHwt8C75kqT4ziLg3efuYZ3n7mm2jKcY2Que1CvArab2Q4zGwO+DVwwLc0FwI0x/D3gPC1051KO4zhORczlMvMTgN0l513AmTOlMbOipAFgGdBTmkjSB4APxNOCpEfnRPHcspxp5aoBalEzuO75phZ116JmqF3dLz6am2riPSgz+wrwFQBJ95vZK6ssadbUou5a1Ayue76pRd21qBlqW/fR3DeXQ3x74LAXlVfFuLJpJGWAVuDAHGpyHMdxaoS5NFD3ARslrZNUB1wEbJ6WZjPw3hi+ELjdzObGqZPjOI5TU8zZEF+cU/oj4FYgDVxnZr+Q9FngfjPbDHwV+Lqk7UAvwYgdia/MleY5phZ116JmcN3zTS3qrkXNsMh0yzssjuM4ThJxTxKO4zhOInED5TiO4ySSxBqoWnSTVIHm90naL+nhePxuNXROR9J1kvbN9H6ZAp+L5XpE0hnzrbGMpiNpPkfSQEldXznfGssh6URJP5P0mKRfSPrjMmkSVd8Vak5cfUtqkHSvpJ9H3Z8pkyaJ7UglupPalqQlPSTp5jLXZl/XZpa4g7Co4v+A9UAd8HPglGlpLge+FMMXAd+pAc3vA/612vVbRvvZwBnAozNcfzPwY0DAWcA9NaD5HODmausso6sTOCOGW4Any3xPElXfFWpOXH3H+muO4SxwD3DWtDSJakdmoTupbclHgG+V+y4cTV0ntQdVi26SKtGcSMzsDsIqypm4APiaBe4G2iR1zo+68lSgOZGY2V4zezCGh4BtBI8qpSSqvivUnDhi/Q3H02w8pq8KS1o7UqnuxCFpFfAW4NoZksy6rpNqoMq5SZr+D3GYmyRgyk1StahEM8A74rDN9yQd7Y7L802lZUsar47DJD+W9NJqi5lOHOI4nfALuZTE1vfzaIYE1ncccnoY2Af81MxmrOuEtCNARboheW3JPwN/CkzOcH3WdZ1UA7VQ+RGw1sxOA37KoV8TzgvPg8AaM3s58Hngh1XWcxiSmoHvAx82s8Fq66mEI2hOZH2b2YSZ/QrBk82rJL2s2poqoQLdiWpLJL0V2GdmD7yQz02qgapFN0lH1GxmB8ysEE+vBV4xT9qOlUo+j0RhZoNTwyRmdguQlbS8yrIAkJQlNPTfNLMflEmSuPo+kuYk1zeAmfUDPwPeNO1S0tqRw5hJdwLbktcCb5f0NGF641xJ35iWZtZ1nVQDVYtuko6oedo8wtsJY/m1wGbg0ri67CxgwMz2VlvU8yHp+KnxbUmvInzXq97wRE1fBbaZ2dUzJEtUfVeiOYn1Lek4SW0xvISwN93j05IlrR2pSHfS2hIz+6SZrTKztYS273Yze8+0ZLOu60R6M7e5c5M0Z1So+UOS3g4UCZrfVzXBJUi6ibAKa7mkLuDThIlZzOxLwC2ElWXbgRHgsuooPUQFmi8E/kBSEcgDF1W74Ym8FrgE2BrnGAD+HFgNia3vSjQnsb47gRsVNk9NAf9uZjcnuR2JVKI7kW3JdI61rt3VkeM4jpNIkjrE5ziO4yxy3EA5juM4icQNlOM4jpNI3EA5juM4icQNlOM4jpNI3EA5c4okK31hT1ImemF+jrfjFzjfGyRdeJT3XhxdyGyVdJekl7/Q+pxDSPqspDdUW4eTPBL5HpSzoMgBL5O0xMzyhJcOE+2FAtgJ/KqZ9Uk6n7Bd9ZlV1lQWSWkzm6hS3pnoU+2YMLOqb83hJBPvQTnzwS0EL8cAvw3cNHVBUpPC3k73Kuwjc0GMXyvpTkkPxuM1Mf4cSVuig8zHJX3zSB6RJZ0Xn7015lUf498cn/GAwv5LNwOY2V1m1hdvv5vgaqjcc4cl/ZWCg9S7Ja0s0X577IXdJml1jL8h5nOXpB1TPbzYg5ja12ePpOtj/HtivTws6cvxxc2pfP9R0s8JDlrLlm+a1pMk/SSW9U5JL4nx/yHp0hj+oKRvxvAWSf8S8340eodA0lWSvi7pfwgvXaYl/b2k+2J5PxjTdUq6o+T+18e0N8TzrZL+pKRepupips/qaUmfid+FrVP6nQVOJXt8+OHH0R7AMHAawb1+A/AwJXsHAX8NvCeG2wh7DTUBjUBDjN9IeBudeO8AwWikgP8FXlcm3xsI3g0aCB6UT47xXwM+XBK/LsbfRPk9bD4GXDtD2Qx4Wwz/HfCpGP4R8N4Y/h3ghyWavht1n0LYnqX0eW3AVoJftU3xOdl47YvApSX5vjOGy5avjNbbgI0xfCbBzQzASoLHitfHuu+I8VuAa2L4bOK+W8BVwAPAknj+gZJy1wP3A+uAjwJ/EePThH2kXkHwzH2wvJV8VjH8NHBFDF8+02fix8I6vAflzDlm9giwltB7umXa5TcCn1BwobOF0EitJrgtukbSVkKjfkrJPfeaWZeZTRIM3trnyf7FwE4zezKe30hocF8C7DCznTH+puk3Svo14P3An83w7DFgai7tgRIdryZs2gbwdeB1Jff80MwmzewxgnGYykvAN4CrLXiEPo/QoN8X6+Y8wmaYABMEx63PV77ScjQDrwG+G5/1ZYI7HcysG7iS4JD0o2ZWusfWTTHNHcBSRf9wwGYLw7UQPr9L43PvIWyfsJHgm/IySVcBp1rYR2oHsF7S5yW9CZjuEf1IZZlyUlta184CxuegnPliM/APhB5Q6R4wAt5hZk+UJo4NWzfwckKPY7TkcqEkPMEcfI8lnUbwEn2+mc3k9HTczKZ8hVWqo1R76dDkVUCXmV1fcu1GM/tkmWeM2uzmnVJAv4XtG8pxKsGx64umxU/3gzZ1niuJE6Fnc+v0h0o6mzC0e4Okq83sawoLTn4D+H3gnYQeZqVM1d2cfOZO8vAelDNfXAd8xsy2Tou/Fbhiah5J0ukxvhXYG3tJlxCGiY6GJ4C1kjbE80uA/4rx6xU24AN419QNcc7oB8AlJb/mZ8NdHHKEeTFw5/MllvQ24A3Ah0qibwMulLQipumQtKbM7TOV7yAW9m7aKem34rMUDcWU5/HzCZsQfkzSupJb3xXTvI7gUX2gTP63EpzEZmPakxXmFdcA3WZ2DcHQn6Gw/UbKzL4PfAo4Y7ZlcRYX/ivEmRfMrAv4XJlLf0nYifMRSSnCCrq3EuZcvh8n8H/C4b/aZ5PvqKTLCMNbGcLQ05fMrCDpcuAnknIxfoorCb28L0a7WTSzV84i2yuA6yV9HNjPkb2Rf4Sw2+i9Mb/NZnalpE8B/xnrZRz4Q2BXJeUrk8fFwL/FZ2aBb0t6HLgGuMzMnpH0UeA6SefGe0YlPRTTz9TTuZYw3PZg/JGxH/hNQk/545LGCfOQl8YyXh/LA3BY73AWZXEWCe7N3Fm0SGo2s+HYsH4BeMrM/qnaupKApC3Ax8zs/mprcRYvPsTnLGZ+L07u/4IwpPjlKutxHKcE70E5juM4icR7UI7jOE4icQPlOI7jJBI3UI7jOE4icQPlOI7jJBI3UI7jOE4i+X9Pf2YXLhhw6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%time gene_mask = utils.select_genes(counts.T, n=3000, threshold=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44808, 3000)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = log_counts.T[:, gene_mask].toarray()\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "x -= x.mean(axis=0)\n",
    "x /= x.std(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min 2s, sys: 1.76 s, total: 2min 3s\n",
      "Wall time: 23.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "U, S, V = np.linalg.svd(x, full_matrices=False)\n",
    "U[:, np.sum(V, axis=1) < 0] *= -1\n",
    "x_reduced = np.dot(U, np.diag(S))\n",
    "x_reduced = x_reduced[:, np.argsort(S)[::-1]][:, :50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44808, 50)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_reduced.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44808,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cell_types.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dict = {\"pca_50\": x_reduced,\n",
    "             \"CellType1\": cell_types.values.astype(str),\n",
    "             \"CellType2\": cluster_ids.values.astype(str)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 732 ms, sys: 16 ms, total: 748 ms\n",
      "Wall time: 755 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with gzip.open(\"data/macosko_2015.pkl.gz\", \"wb\") as f:\n",
    "    pickle.dump(data_dict, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
